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1 Final Report 

The primary Phase I research and development objectives were to: 

1. Develop curvature-driven, flow-based multiscale image processing methods. We pro¬ 
posed to develop image enhancement and segmentation methods using nonlinear par¬ 
tial differential equation based models and algorithms, level-set techniques, and active 
contour models. 

2. Integrate the proposed methods into an open and extensible software environment. We 
proposed to integrate the developed methods into our open and extensible software en¬ 
vironment featuring S-Plus and Java-based client tools for processing and displaying 
images. 

3. Demonstrate the proposed methods on medical and tactical data. We proposed to apply 
combinations of the methods developed to specific applications in medical and tactical 
imaging. 

The primary Phase I results are: 

1. MathSoft has focused on the development of level-set and variational level-set segmen¬ 
tation techniques based on [1-5,16]. 

2. MathSoft has implemented these algorithms in C and linked them to the S-Plus 
language to enable scripting and simulation of complex applications. Our implemen¬ 
tations can process both two- and three-dimensional image data. 

3. MathSoft has applied these algorithms to the segmentation of the surface of the brain 
from a two-dimensional MRI dataset and the segmentation of the prostate from two- 
dimensional ultrasound images. 

4. Level Set Systems, the Phase I subcontractor, has focused on the development of an 
active contours technique for curve evolution [6] that makes use of a Mumford-Shah [7] 
like fitting term. 

A detailed discussion of the research conducted and our results is given in this section. 
We begin with a general overview of level-set methods in Section 2 and its application 
to image segmentation in Section 3. Implementation of level-set methods requires careful 
attention to accuracy and robustness as well as computational speed. In Section 4, we 
discuss the implementation issues we have encountered and solutions to them. Preliminary 
results are given in Section 5. Simple segmentation methods often fail in images where there 
are indistinct or missing boundaries. In Section 6, we discuss two incremental improvements 
we used to address this problem. In Sections 7 and 8, two new models using region-based 
features are introduced which provide further robustness to weak or blurred edges. 
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2 Introduction to Level-Set Methods 


The level-set method developed in [8] has proven to be phenomenally successful in repre¬ 
senting and following the motion of interfaces or surfaces, not only in image processing, 
computer vision and graphics, but also in fluid dynamics, materials science, combustion, 
and elsewhere. This approach represents a surface as a level-set of a smooth function, and 
moves the surface solely by evolving the representing function via a nonlinear partial differ¬ 
ential equation (PDE). The velocity of the motion often involves curvature. This “implicit” 
representation of a surface turns out to have major advantages over traditional methods, 
especially in its ability to handle topological changes such as merging and pinching off. This 
formulation of the surface motion and the associated numerical methods were abstracted 
from techniques developed originally for combustion and fluid dynamics simulations. 

We now describe briefly the basic level-set approach. The idea is to deflne a smooth level- 
set function (^(x, t) that represents the evolving surface r(t) as the set where ip{x, t) = 0, with 
(p{x, t) < 0 representing the interior bounded by the surface and (p{x, t) > 0 representing 
the exterior. The evolving surface is captured for all time with r(T) = {^{x,t = r) =0}. 
This deceptively trivial statement is of great significance because topological changes such 
as merging and breaking are well defined, and are performed automatically and implicitly. 

The surface motion is modeled by convecting (p with the desired surface velocity field v 
(and arbitrary elsewhere) via the equation 

^+vS/p = 0. ( 1 ) 

Actually, only the normal component of v is needed Vn = v ■ |^, so that Eq. (1) becomes 

^+u„|V^| = 0. (2) 

In Phase I research, we focused on variational level-set methods and their application to 
the segmentation of the surface of the brain from a three-dimensional MRI dataset. Many 
physical systems evolve to minimize an energy functional, and the variational formulation 
leads to a very stable discretization. As mentioned above, these ideas have been used 
successfully in image processing [9], computer vision [1], and related fields. The energy 
functional, which involves surface energy and bulk energy, can be represented simply and 
completely by level-set functions. 

Among the advantages of variational level-set methods are: (I) A numerically stable and 
robust algorithm is guaranteed; (II) a natural physical interpretation can be incorporated 
into the formulation; (III) constraints such a fixed volume, boundary conditions (based 
on a penalty method), and normals, curvatures, and other geometric quantities may be 
prescribed; (IV) shapes may be deformed dynamically by coupling to external physical 
laws [10]. 
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3 Level Sets for Image Segmentation 

The level-set method has been applied to the problem of image segmentation [2]. The basic 
idea is to start from an initial seed and propagate a contour or front r{t) until the contour 
stops at the shape boundary. This is accomplished by defining a speed function based 
on image gradient features. Note that the propagation of the contour is not done directly. 
Rather, it is done through the evolution of a level-set function (p. The advantages of this 
method are firstly that arbitrarily complex shapes (e. g. corners and protrusions) can be 
recovered, and secondly that more than one shape can be isolated at a time as topological 
changes are handled implicitly. 

In [2], the speed function has the form 

Vn = ki {±1 — ck) (3) 

where k is the local curvature of the level-set. The parameter e controls the smoothness of 
the propagating front, and hence its robustness to noise. The sign of the first term depends 
on whether we wish to propagate the interface outwards (positive) or inwards (negative). 
The quantity k[, known as the image speed term, depends on the image gradient 

" l + \VG,*I{x,y)\ 

where {G^ * I) denotes the convolution of the image with a Gaussian smoothing filter, kj is 
close to one in regions of low image gradient and close to zero in regions of high gradient. 
Hence, ki forces the propagation to stop near the shape boundary where we expect there to 
be a large change in the image intensity. 

Evolution of the level-set function p is governed by Eq. (2). We discuss the implemen¬ 
tation of this evolution process in the next section. 


4 Implementation of Level-Set Methods 

4.1 Numerical Approximation 

To prevent smoothing of discontinuities, such as corners, care is needed in the numerical 
approximation of the evolution equation, Eq. (2). Several entropy-satisfying schemes are 
outlined in [5]. These make use of upwind schemes where derivatives are calculated using 
only neighboring pixels from the direction of propagation. 

4.2 Signed Distances Calculation 

The level-set approach requires an initial function ip{x, t = 0) with the property that the 
zero level set {p = 0} corresponds with the initial front r(0). A straightforward and 
computationally expensive technique is to calculate the signed distance function from each 
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pixel, or image grid point, to the initial front. Note that the initial front r(0) is defined on 
a continuous domain not just at grid points. As some level-set evolution schemes require 
periodic re-initialization of the level-set function to the signed distance function, a fast 
and accurate approximation method is needed. In our research, we investigated two signed 
distance approximation schemes: the distance transform operator [11] and the Fast Marching 
method [4]. 

Distance Transform Operator The distance transform operator [11], has the advantage 
that it can quickly approximate the distance from each image grid point to the zero level set 
since it requires only two passes through the image. We found that the major disadvantage 
of this method is that it assumes that points on the continuous zero level set are located 
only at regularly spaced lattice grid points. Approximating a point on the zero level set by 
the closet grid point results in significant loss of accuracy, particularly when re-initialization 
is done regularly. 

Additionally, to use the distance transform for re-initialization, we have to explicitly 
extract the zero level set from the continuous function cp. One method is to construct a 
polygonal approximation from the grid cells that contain a zero crossing [2]. This requires 
a contour tracing procedure to link neighboring cells to form a closed polygon. Contour 
tracing is a difficult problem due to effects of noise, sharp protrusions and topology changes. 
We implemented a tracing procedure that effectively performs a depth-first search of the 
pixel/cell space. Although we have attempted to cover most contingencies, robustness is not 
guaranteed. Protrusions may be smoothed, and in some cases smaller contours are retained 
at the expense of larger contours being found. 

Due to errors introduced in both the distance transform and contour extraction, we 
do not recommend this scheme for use in level-set implementations. In our experiments we 
found that these errors caused the zero level set to be perturbed away from shape boundaries. 

Fast Marching The distance calculation from the zero level set can be viewed in terms 
of front propagation [4]. This simplifies to solving the following Eikonal equation 

|vr| = 1 

where T{x, y) is the arrival time of the front at image grid point (a:, y). The Fast Marching 
method [5] is a fast technique for solving Eikong,! equations. Fast Marching requires only 
one pass through the grid points since it takes advantage of the causality relationship of the 
front propagation. 

The difficulty lies in the initialization stage where we need to provide approximate dis¬ 
tances at a set of trial points to begin the Fast Marching. We use an initialization method 
that tags as trial all image grid points inside the front with at least one neighbor outside 
the front [4]. The approximate distance is calculated by linearly interpolating the level-set 
function value between the grid point and its neighbors to find the distance to the zero set. 

The major advantages of this method are that it has fast computation and that the 
constructed signed distance has the same zero level set as the given level-set function. We 
have implemented this method in both in two and three dimensions, and have used it in all 
our segmentation algorithms. 
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4.3 Image Speed Extension 

The image speed term ki in Eq. (3) has meaning only on the front T{t), as it was designed to 
force the propagation to stop near the shape boundary. However, to solve for (f in Eq. (2), 
ki must be defined for all grid points. ki can be extended as follows [2]; for every point P 
not on the front T{t), ki{P) is assigned the value of ki{Q) where Q is the closest point to P 
that lies on the front. This extension method is computationally expensive and forms the 
major bottleneck when executing the algorithm. 

The extension oi kj outward from the zero level set is a problem similar to that of the 
signed distance calculation from the zero level set. In both problems, we need to locate at 
every grid point the closest point on the zero level set. 

The Fast Marching method we used to calculate the signed distance can also be used to 
simultaneously extend kj or any other quantity [4]. During the fast marching, a grid point 
is assigned a /c/ value that is the weighted average of the kj values at the grid points used 
to compute the distance. 

In the initialization phase, if ki is defined on a finer grid than the calculation grid, sub¬ 
grid resolution can be taken into account to define the kj values at the trial calculation grid 
points. These values are then extended smoothly out by the Fast Marching process. Another 
advantage of this method is that ki is extended such that the signed distance function is 
maintained as the level-set function evolves, preventing the level sets from bunching up or 
spreading out. 

Our implementation of the Fast Marching method can handle seamlessly both two- and 
three-dimensional datasets and can extend any number of quantities at the same time as 
it calculates the signed distance. We have used this implementation in all of our level-set 
segmentation algorithm prototypes. 

4.4 Level-Set Evolution Algorithms 

The flowchart of the basic level-set evolution algorithm is shown in Figure 1. In our imple¬ 
mentation, Fast Marching is used in step (1) to calculated the signed distance function and 
to extend the image speed term kj. ^ 

A more efficient algorithm can be implemented where the level-set function is updated 
only at a small set of points within a band of width 5 around the zero level set of [2]. 
The width 5 of the narrow band determines the number of points that need to be processed. 
The zero level set that lies inside the narrow band moves until it collides with the boundary 
of the narrow band. Upon collision, the level-set function (p is re-initialized by treating the 
zero level set of as the initial contour r(0). 

This narrow-band algorithm requires the selection of a parameter Z, which is the number 
of time steps required to move the front by a distance of roughly 5/2. The choice of I requires 
some experimentation. Alternatively, we can re-initialize p when we detect a possible zero¬ 
passing through a grid point outside of the narrow-band. A flowchart for the narrow-band 
scheme is shown in Figure 2. 
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Figure 1; Flowchart of the basic level-set evolution algorithm. 


A third alternative is to re-initialize the level-set function at every iteration. Thus, the 
zero level set will always lie in the center of the narrow-band, and hence collision detection 
with the narrow-band boundary is not required. Re-initialization at every step is achieved at 
no extra computational expense as we are already required to extend kj at every step. Using 
the Fast Marching technique both are done simultaneously. This scheme is diagrammed in 
Figure 3. 

Identification of the narrow-band is another important implementation issue. Options in¬ 
clude using morphological operations or moving along the zero level set and collecting points 
within a distance of (5/2. Another alternative is to'modify the Fast Marching re-initialization 
method to collect the coordinates of each point being processed, and to terminate the march¬ 
ing when the distance has reached a value greater than 5/2. We have used the latter method 
in our algorithm prototypes. 

4.5 Algorithm Framework 

We have developed an extensible modular framework for level-set implementations. The 
framework can seamlessly handle two- and three-dimensional data. New algorithms can be 
rapidly prototyped by building new modules to extend the existing set of modules. For 
efficiency and portability the software was developed in C. 
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(a) Narrow-band algorithm with updating. (b) Update algorithm. 

Figure 2: Flowchart of the narrow-band algorithm, (a) The overall algorithm, (b) The 
level-set function update algorithm. 

5 Preliminary Results 

5.1 Application to Synthetic Data: Single Circle 

This simple test image, Figure 4(a), is of a uniforin intensity circle on a black background. 
The initial contour. Figure 4(b), is small square located within the circle, below and to the 
right of the center of the circle. The zero set expands uniformly until it reaches the boundary 
of the circle, Figure 4(c). The contour near the boundary stops expanding while the rest of 
the contour continues to grow until it reaches the rest of the circle boundary, Figure 4(d). 

5.2 Application to Synthetic Data: Two Circles 

For this test the image contains two circles of different radii. Figure 5(a). The initial contour 
is a square that encloses both circles. Figure 5(b). In this example, the zero set contracts 
uniformly until it reaches the boundary of the circles. The zero set continues to contract 


8 













Figure 3; Flowchart of the narrow-band algorithm with re-initialization at every iteration. 

around both circles until the contour meets up with itself between the circles, Figure 5(c). 
At the point, the zero set changes topology and splits into two contours. The contours 
continue to contract until they reach the boundary of each of the circles. Figure 5(d). 

5.3 Application to MRI Medical Imagery: Brain 

This test image is a cross-section of an MRI brain scan. Figure 6(a). The initial contour 
is a thin rectangle located within the gray matter, Figure 6(b). Figure 6(c) shows the 
results obtained using the distance transform operator method for re-initialization. Note 
that the sharp indentation towards the bottom of the shape has been smoothed. In contrast. 
Figure 6(d) shows the results obtained using the Fast Marching re-initialization, which 
captures accurately the indentation. 
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Figure 4: (a) Image of a uniform intensity circle, (b) The initial contour 7 ( 0 ) is a small square located within the 
circle, below and to the right of the center of the circle, (c) The zero set at iteration 10. The boundary has expanded 
uniformly to the boundary of the circle, (d) The zero set, at iteration 18. The contour already at the boundary of the 
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Figure 5: (a) Image of two circles of different radii, (b) The initial contour 7 ( 0 ) is a square that encloses both circles, 
(c) The zero set at iteration 15. The zero set contracts uniformly until it reaches the boundary of the circles. It 
continues to contract around both circles until the contour meets up with itself between the circles, (d) The zero set 
at iteration 21 . The zero set changes topology and splits into two contours. The contours continue to contract until 
they reach the boundary of each of the circles. 













(a) Original image. 



(c) Result when using polygonal 
approximation and distance trans¬ 
form. 



(b) Initial contour. 



(d) Result when using the Fast 
Marching re-initialization. 


Figure 6 : (a) MRI image of a brain, (b) The initial contour 7 ( 0 ) is a square inside the 
gray matter, (c) The zero set obtained using the polygonal approximation and distance 
transform. Note that perturbation of the zero set at each re-initialization has caused the 
indentation at the bottom of the image to be smoothed, (d) The zero set obtained using 
the Fast Marching re-initialization. 


5.4 Application to Three-Dimensional Synthetic Data: Two Spheres 

In this example, we wish to segment two spheres in a three-dimensional data set with 
dimensions 60 x 60 x 60. The origin, ( 1 , 1 , 1 ), is placed on the lower left-hand corner. 
Each of the spheres has a radius of 5 pixels and are located at positions (20,20,20) and 
(40,40,40). Figure 7 shows the zero set at different times during the propagation process. 
The final images were constructed by building a mesh over the surface using the method 
marching cubes provided in a public domain software called isosurf. We have rendered the 
mesh using another public domain program called Geomview. These program are available 
at: 

http: //svr-www. eng. cam. ac. uk/*'gmt 11 /software/isosurf/isosurf. html 
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http;//www.geom.umn.edu/software/geomview/ 

The initial contour, Figure 7(a), is a sphere of radius 26 pixels located at (30,30,30). As the 
propagation process continues, the zero set contracts uniformly until it reaches the boundary 
of the spheres. Figure 7(b). The zero set continues to contract around both spheres until 
the surface meets itself between the spheres. Figure 7(e). At this point, the zero set changes 
topology and splits to form three spheres. Figure 7(f). Further contraction of the two outer 
spheres is stopped by image speed ki, while the center sphere continues to contract to a 
point. Figure 7(h). 
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zero set at iteration 25. The parts of the zero set around the spheres boundary stops contracting while the rest of the 
surface continues to shrink, (d) The zero set at iteration 30. (e) The zero set at iteration 35. The zero set continues 
to contract around the spheres until it meets itself, (f) The zero set at iteration 40. The zero set splits into three 
surfaces, (g) The zero set at iteration 45. Contraction of the two outer surfaces stops due to the boundaries of the 
spheres, while the center surface continues to contract, (h) The zero set at iteration 54. The converged solution. 













6 Improvements 

Simple segmentation methods may fail in images where there are indistinct or missing bound¬ 
aries. In this section, we explore two incremental improvements to address this problem: 
the incorporation of a priori edge information, and the use of geodesic active contours. 

6.1 Incorporating A Priori Edge Information 

An example of indistinct boundaries is shown in the ultrasound image of a prostate in 
Figure 9(a). By calculating A:/, the image speed term, based solely on image gradient infor¬ 
mation, Figure 9(b), it is difficult to distinguish between the prostate and the background. 
For this image, the level-set algorithm fails to stop at the shape boundary. 

To remedy this, we have used a priori edge information obtained from edge detection 
preprocessing (see the flowchart in Figure 8). The edge image in Figure 9(c) is obtained 
by using anisotropic diffusion smoothing [12], followed by Canny edge detection [13]. Fig¬ 
ure 9(d) is the image speed term {kj) determined using the thresholded square distance 
to the closest edge pixel. Using this ki, a good outline of the prostate can be obtained 
(Figure 10). 

Note that while the original image has dimensions 512 x 512, the Fast Marching extension 
method allows us to obtain accurate results using sparser calculation grids. The results 
shown in Figure 10 were obtained using a calculation grid of 256 x 256 and 128 x 128 points. 
This reduction in grid size provides a large saving in computation time. 
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Original Implementation - using no a priori information 


Using a priori edge information 



Figure 8: Flow charts of the level-set segmentation process with and without a priori infor¬ 
mation. 
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6.2 Geodesic Active Contours 


A partial solution to the indistinct boundary problem can be found in [14,15], where the 
propagation speed function is given as; 


Vn^ ■ 



V kj -V (p 
IIV^II 


— kiK. 


(5) 


The second term in Eq. (5) acts as smoothing term to reduce the length of the contour. The 
first term acts like a doublet which attracts the contour to the potential wells in kj. This 
technique is known as a geodesic active contour in [14], and conformal length shortening 
in [15]. 

An additional constant inflation (deflation) term is also suggested to speed up conver¬ 


gence [14,15]: 


Vn 


V kj -V (f 

liv^ll 


+ ki{c - k). 


( 6 ) 


Note that the image speed terms in Eq. ( 6 ) have meaning only on the front r(t). Hence 
at each iteration of the algorithm we need to extend all components of Vkj and well as kj. 
This additional extension requires only a small amount of extra work in our Fast Marching 
re-initialization method. 

The major differences between using the speed functions in Eq. (3) and Eq. (5) are 
illustrated in the following simple example that uses a synthetic test image consisting of a 
uniform intensity circle on a black background. In this example. Figure 11, we use an initial 
contour 7 ( 0 ) of a circle that is the same size as the one we wish to recover, but we shift it 
to the upper left-hand corner by a few pixels. 

The results obtained using the original speed function, Eq. (3), are shown in Figure 12, 
which illustrates the behavior of the old speed function once the shape boundary has been 
passed. Using a constant inflation term in Figure 12(b), the contour expands to meet the 
circle boundary except in the top-right region where it leaks into the background. Using a 
constant deflation term in Figure 12(c), the contour contracts to meet the circle boundary 
except in the bottom-left region where it shrinksJnside the shape. 

For the same initial contour. Figure 13 shows the results obtained using the new speed 
function, Eq. (5), which allows movement both inwards and outwards to try follow the bound¬ 
ary of the circle. Note that the bottom-right boundary is slow to converge. Figures 13(b) 
and 13(c). This is partially due to the curvature term acting to pull the contour away from 
the boundary. In this example, convergence does not occur until the 57th iteration. The 
rate of convergence can be improved by adding a constant inflation component to the speed 
function as in Eq. ( 6 ), where c> 0. With c set to 0.05 for the same test, convergence was 
reached in 37 iterations, a significant improvement when compared to 57 iterations when 
c = 0 . 
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(a) Test image. 


(b) Edge potential kj. 




(c) Column derivative dkjfdx. 


(d) Row derivative dkijdy. 


Figure 11: (a) A simple test image consisting of a uniform intensity circle on a black back¬ 
ground. (b) The image of the edge potential kj calculated using Eq. (4). The grayscale 
values range from black (0) to white (1). Images of the (c) column and (d) row derivatives 
of kj. The grayscale values range from black (-0.5) to white (0.5). 
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Figure 13: Segmentation using speed function Eq. (5). (a) The initial contour 7 ( 0 ) is a circle that is the same size as 
the circle we wish to recover but offset by few pixels, (b) The zero set at iteration 20. The contour shifts to meet the 
boundary, except in the bottom-right region, (c) The zero set at iteration 40. The contour converges slowly near the 
bottom-right region, (d) The zero set at iteration 57. The converged solution. 












7 Segmentation Using Region Statistics 

In Section 5, we observed that edge-based segmentation methods fail when there is insuf¬ 
ficient intensity change at the object boundaries. Two incremental improvements of incor¬ 
porating a priori edge information and using geodesic active contours were discussed in 
the last section. Another solution is to use region-based features. The use of region-based 
features avoids the need to calculate image gradients that are extremely sensitive to noise, 
and provides greater robustness to the effects of weak or blurred edges. 

For example, the region-based method in [16] assumes that an image consists of a finite 
number of regions, where each region is delineated by a predetermined set of features or 
statistics (e. g. means, variances, and textures). In [16], an initial curve r(0) is evolved to 
maximize the difference between the value of the feature(s) inside the curve from the value 
on the outside. The curve evolution is implemented in a level-set framework, where r(t) is 
embedded as the zero level set of the level-set function ip. 


7.1 Bimodal Images Segmentation 

Suppose we wish to segment a bimodal image, i. e. an image that consists of two regions, 
possibly multiply connected, delineated by a particular feature or statistic. By multiply 
connected, we mean that each region may consist of many disconnected parts; for example, 
this page consists of two regions, one white and one black, where each region has many parts 
that are disconnected, such as letters. The approach of [16] is to minimize the following 
energy functional: 

E(r) = -5(»-v)^ (7) 

where u is the value of the feature for the region inside the curve F and v is the value for 
the region outside F. At each iteration, F is evolved in a steepest descent manner such that: 


dt 


(u — v) (Vu — Vu). 


If u and V represent the mean intensities, then we have: 


dt 


, . (I — u I — v\ 


N 


( 8 ) 

(9) 


where I is the image intensity. An, Ay are the area of the regions inside and outside F and 
N denotes the outward unit normal of F. This evolution can be translated to a level-set 
evolution such that: 


dip 

dt 



I — V 

Ay 


\Vip\ = 0. 


( 10 ) 


The use of a level-set framework allows complex shapes to be segmented and topological 
changes such as merging and splitting to be handled implicitly. Since this approach uses 
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region-based statistics, the initial curve r(0) can be placed anywhere on the image domain 
(e. g. the initial curve may overlap object boundaries). 

To add robustness to noise, a curvature-based term can added to second term of Eq. (10): 

+ = ( 11 ) 

Note that this approach does not require the extension of the image speed in contrast 
to other methods described in this report. However, in our experiments we have found that 
the level sets may bunch up making the curvature calculation noisy. Therefore, regular re¬ 
initialization of the level-set function y? to the signed distance function (from the zero level 
set) is recommended. 

As described in [16], more than one feature can be used to segment a bimodal image. 
The approach is then to maximally separate the two feature vectors using some appropriate 
distance measure. 

Extension to trimodal images is also described in [16]. To segment an image into three 
regions/classes at least two features are needed. Two (coupled) level sets are evolved to 
maximize the area of the triangle formed by the feature vectors. In general, to segment an 
image into N classes, at least N — 1 features and N — 1 level sets are require. 

7.2 Implementation Issues 

When dealing with an image, updating the level set at every image grid point can be 
computationally intensive. To improve efficiency, the narrow-band scheme described in 
Section 4 can also be applied here. In this scheme, the level-set function is updated only at 
a small set of points within a band of width S around the zero level set of <p". 

In our narrow-band implementation, cp’^ is re-initialized to the signed distance function 
at every iteration. This strategy means that the zero level set will always remain in the 
middle of the narrow-band. Therefore collision detection with the narrow-band boundary is 
not required. 

In the narrow-band scheme, the choice of the time step At is important. If At is too 
small, the algorithm will be slow to converge; if At is too large, the zero level set may move 
out of the narrow-band. A good choice of At is the one that keeps things moving if a large 
proportion of the zero level set is in a region of approximately constant intensity. In our 
experiment, we found a satisfactory choice of At is: 

^ ^ _ min (Au, A„) 

. \2 ■ 

\U — v) 

Note that in images where the region we wish to segment is multiply connected, the 
narrow-band scheme may miss regions which are never within 5 of the zero level set. 
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7.3 Application to Microscopy Imagery: Red Blood Cells 

Figure 14(a) shows a microscopy image of red blood cells. The number and size of red blood 
cells (RBC) in a sample can be used to diagnose anemia and other diseases. In this example, 
we apply the region statistics method [16] to segment the RBCs from the background. 

The initial curve is a circle located in the top-right region of the image, Figure 14(a). 
The curve evolves to separate maximally the mean intensity inside the curve and the mean 
intensity outside the curve. Note that the evolving curve can move both in and out. For 
example, consider the RBCs near the initial curve. In Figure 14(b), it can be seen that the 
zero level set has move inwards to locate the RBCs inside the initial curve as well as locating 
RBCs outside the initial curve. 

Figure 14(c) shows the zero level set after 800 iterations. This example illustrates the 
major advantage of using level sets: topological changes are handled automatically allowing 
all RBCs to be segmented. The level-set framework also makes it simple to extract the 
RBCs for post-processing. The RBCs are regions where cp > 0. These are shown as dark 
grey regions in Figure 14(d). The dark gray value is the steady state mean for regions 
outside r, while the light gray is the steady state mean for regions inside F. 

Individual RBCs can be separated using connected component labeling [11]. In Fig¬ 
ure 15(a), each individual RBCs is assigned a different gray value. Note that we have also 
eliminated all RBCs that are not wholly within the image. 

From Figure 15(a), it can be seen that except for two cases, the segmentation method 
was able to separate RBCs which were very close together in the original image. Using 
the image in Figure 15(a), simple statistics can be calculated. Figure 16. The image in 
Figure 15(b) shows the smallest and largest cell detected. The small “cell” is actually a 
platelet and the large “cell” is really two cells close together. 

7.4 Application to Prostate Boundary Segmentation 

Measuring the size of the prostate and identifying accurately the prostate boundary is an 
important part in the diagnosis and treatment of prostate cancer. Three-dimensional trans- 
rectal ultrasound (TRUS) imaging is routinely used for prostate examination. Ultrasound 
images are typically low contrast and images of soft tissues (e. g. the prostate) suffer from 
blurred boundaries. In this section, we apply the region statistics method of [16] to the 
problem of prostate boundary segmentation from TRUS images. 

Results using two-dimensional segmentation are shown in Figures 17 and 18. Figure 17(a) 
shows one of the center slices of a TRUS image set. The initial curve is a large circle centered 
around the center of the prostate. The curve evolves in a manner that tries to maximize 
the difference between the mean intensity inside the curve and the mean outside the curve. 
Figure 17(b). 

Figure 17(d) shows the zero level set after 200 iterations. It can be seen than the general 
shape of the prostate has been captured. The ragged edges on the bottom-left are due to 
the blurred boundaries in the original image. These edges can be smoothed out in post- 
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processing. However, too much smoothing may lead to the underestimation of the prostate 
size. The size of the prostate can be estimated by counting the number of pixels where 
tp <0. For this image, the prostate is approximately 37,298 pixels in area. 
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(a) Initial contour. 



(c) Zero set, iteration 800. 


B'as, 

""Sit* 



. 

P'Qik 

(b) Zero set, iteration 400. 



(d) Segmented red blood cells. 


Figure 14: Red blood cells segmentation from a 169 x 169 microscopy image, (a) The initial 
contour r(0) is a circle near the top-right corner of the image, (b) The zero level set at 
iteration 25. The contour has evolved to maximally separate the mean intensity inside the 
contour from the mean intensity on the outside^ Note that the contour has split several 
times to segment red blood cells individually, (c) The zero level set at iteration 800. All 
the red blood cells have been isolated, (d) The reconstructed image. The dark gray region 
represents the region inside the contour of image (c). The dark gray value is the mean 
intensity of all pixels inside the contour, while the light gray value is the mean intensity of 
all pixels outside of the contour. 
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(a) Connected component labeling 


(b) Smallest/largest segmented cell. 


Figure 15: Connected component analysis of Figure 14(d). (a) The result from connected 
component labeling. Each disjoint region is represented by a different gray value. Note that 
all regions not wholly within the image has been removed, (b) The smallest (white) and 
largest (gray) segmented regions. The small region is actually a platelet. The large region 
is really two red blood cells that are close together in the original image. 



Number of Cells = 33 
Average Size = 246.8 pixels 
Median Size = 240 pixels 
Smallest Cell = 13 pixels 
Largest Cell = 520 pixels 


Figure 16: Flowchart of the red blood cell counting process. 










(a) Initial contour. 


(b) Zero set, iteration 25. 




(c) Zero set, iteration 50. 


(d) Zero set, iteration 200. 


Figure 17; Prostate segmentation from a 280 x 280 ultrasound image. This is the sixth slice 
of an eight slice trans-rectal ultrasound image set. (a) The initial contour is a large circle 
centered near the center of the prostate, (b) The zero set at iteration 25. The contour has 
evolved to maximally separate the mean intensity inside the contour from the mean intensity 
outside of the contour, (c) The zero set at iteration 50. (d) The zero set at iteration 200. 
The ragged edges at the bottom-left are due to blurred edges in the original image. 
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Segmentation results for other images in the TRUS set are shown in Figure 18: 

• Slice 3 (Figure 18(e)) - interference from an artifact has cause the segmentation to 
also capture the dark region at the top. 

• Slice 4 (Figure 18(f)) - most of the prostate has been captured. There is slight 
underestimation on the top-right. This is due to the region being lighter in color. 

• Slice 5 (Figure 18(g)) - again most of the prostate has been captured with underesti¬ 
mation along the top boundary. The dark fringes on the right have caused “tails” in 
the segmentation results. 

• Slice 7 (Figure 18(h)) - the prostate is barely visible in the original image. The 
segmentation algorithm was able to capture most of the prostate with some underes¬ 
timation at the top. 

Theoretically, all images in a TRUS set can be used together to reconstruct the prostate 
in three dimensions. In practice, this is a difficult problem as there are typically only eight 
slices in/a TRUS image set, therefore resolution in the third dimension is poor. To address 
this problem, intermediate slices can be interpolated from the original images. However, the 
images needs to be downsampled to reduce the huge discrepancy in the grid spacing between 
the dimensions. 

In our experiment, we have linearly interpolated ten intermediate slice between each pair 
of original images. The images were also downsampled to give a dataset of size 70 x 70 x 45. 
Note that we have use only slices 3 to 7 of the TRUS set, as the prostate is barely visible in 
the start and end slices. We then apply the region-based segmentation method of [16]. 

The initial surface r(0) is a partial sphere. Figure 19(a). The surface evolves to maximally 
separate the mean intensity inside from the mean intensity outside, Figures 19(b)-19(d). The 
object at the top is due the dark semi-circle which appears in all the images. 

From Figure 19(d) some of the major characteristics of a prostate are visible. The 
prostate is approximately as long as it is wide and roughly conical in shape. The bumps of 
the two lateral lobes at the top of the prostate can be seen in the reconstruction. 

Figure 20 shows the contours of the surface in Figure 19(d) overlaid on the original 
downsampled slices. From these images it can be seen most of the problem with the 3D 
segmentation are due to blurry nature of slices 3 and 7. Although far from ideal these results 
shows promise and we propose to investigate the use of this region-based segmentation 
method as part of Phase II work. 
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(a) Initial contour. 


(b) Zero set, iteration 5. 



(c) Zero set, iteration 20. 


(d) Zero set, iteration 100. 


Figure 19: Three-dimensional segmentation of the prostate from a set of trans-rectal ul¬ 
trasound images, (a) The initial surface is a partial sphere centered near the center of the 
prostate, (b) The zero set at iteration 5. The surface evolves to separate the mean intensity 
inside the surface from the mean intensity on the outside, (c) The zero set at iteration 20. 
(d) The zero set at iteration 100. The object at the top is due to the dark semi-circle which 
appears at the bottom of all images. The ragged edges at the front are due to the blurriness 
of the image in slice 7, Figure 18(d). 







(a) Slice 3. 


(b) Slice 4. 


(c) Slice 5. 



(d) Slice 6. (e) Slice 7. 


Figure 20: Three-dimensional segmentation of the prostate from a set of trans-rectal ultra¬ 
sound images. Each image shows the contour of surface in Figure 19(d) onto each of the 
original downsampled slices. 
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8 Segmentation Using Active Contours without Edges 
- LSS. Inc. 

Recently, a new model [17] for active contours has been proposed that is based upon using a 
Mumford-Shah [7] like functional. In contrast to classical active contours models or snakes, 
the stopping term in this method has no dependence upon image gradient or edge strength 
V/. Rather, the stopping term is related to the segmentation of the image. This new method 
is advantageous in that objects can be detected that are not necessarily defined by sharp 
edges. Furthermore, this method lends itself readily to a level-set formulation, wherein the 
evolving curve F is equated to the zero level set of a function ff. 

For a given image /, defined on a domain Q and formed by two regions of piecewise 
constant intensities, the main idea behind this algorithm is to minimize an energy based 
segmentation. The key idea is the formulation of a fitting energy 



where the constants Ci and are averages of the image I inside and outside the boundary 
curve r,' respectively. Clearly, the above term is minimized when the evolving curve F 
perfectly delineates the regions formed by the two intensities. Additional regularizing terms, 
dependent upon the length and area within F, may be added to the fitting energy. All terms 
may, in fact, be rewritten in terms of the level-set function (p and the resulting general fitting 
energy is: 

■F((^,Ci,C2) = 

H / 5((p)\V(p\dx + 1 ' / H{(p)dx 
J ^ */ nt 

-k [ \I-c^\'^H{p)dx+ [ \I-C 2 \\l-H{p))dx 
Jn Jn 

where H and 5 are the Heaviside and Dirac functions and p, u are constant fixed parameters 
for the length and area terms, respectively. The associated Euler-Lagrange equations for p 
can be derived by minimizing F with respect to (p and parameterizing the descent direction 
by an artificial time. 

This “active contours without edges” model has the advantage of detecting interior con¬ 
tours automatically. An example of this is shown in Figure 21. The original image for these 
figures is shown in Figure 21(a). The three distinct objects within this image are successfully 
captured by the active contours without edges model. In Figure 22, it can be seen that the 
new model works well on the same image with the addition of noise. 

The active contours without edges model can be easily extended to three dimensions. 
Some examples of this are shown in Figures 23 and 24. In Figure 23, the algorithm has been 
applied to an image that contains two balls. In Figure 24, the algorithm has been applied to 
an image a cube with a semi-spherical region cut out from it. In both examples, the initial 
zero level set of is a large sphere. 
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Figure 25 shows the segmentation of the prostate from a three-dimensional transrectal 
ultrasound image (TRUS). A two-dimension cross-sectional snapshot, Figure 25(a), shows 
the prostate. The original image was cropped to remove extraneous text. 

Figure 26 and Figure 27 show the segmentation of the brain fro a three-dimensional 
image. The grid size use for Figure 26 is 80 x 80 x 7, while the grid size use for Figure 27 
is 29 X 29 X 19. 

Again, the algorithm performs well, as many of the grooved details are captured by the 
(blue) evolving surface. 
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(a) Original image. (b) Iteration 600. (c) Iteration 1000. (d) Iteration 1900. 

Figure 22; Results from the active contours without edges model, (a) The original noisy image contains three 

The initial contour is a circle that encompasses all three objects. Segmentation results after (b) 600, (c) 1000, and (d) 

1900 iterations. 
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The contour after 3000 iterations. 























37 


Figure 27: Results from the active contours without edges model on a three-dimensional brain image using a c 
grid of 29 X 29 X 19. (a) The initial contour. Segmentation results after (b) 60, (c) 100, and (d) 200 iterations. 
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